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1. Introduction 

Systems of particles interacting through long-range forces are usually described over 
a certain time scale by Vlasov equations. This situation is encountered in various 
fields of physics: plasma physics, self-gravitating systems, wave-particles interactions 
for instance. One may add here two dimensional fluid dynamics, since the 2D Euler 
equation shares many properties with the Vlasov equation. 

A Vlasov equation usually admits a continuous infinity of stationary states. 
Investigating the stability of these states is a natural question. In a celebrated paper 
PQ, Landau considered stationary states of a plasma which are homogeneous in space, 
and addressed the issue of the asymptotic behavior of a small perturbation through 
a Laplace transform analysis of the Vlasov equation linearized around the stationary 
state. This was the starting point of an extremely abundant research on Landau 
damping in plasma physics, and more generally on the fate of perturbations around 
stationary states of the Vlasov equation, usually homogeneous in space. We focus now 
on the known mathematically rigorous results, all of them obtained in the context of 
homogeneous stationary solutions of the Vlasov equation. A rigorous linear treatment 
"a la Landau" is provided for instance in [2j [3] : under very strong regularity hypothesis 
for the stationary state and the perturbation (both should be analytic functions of the 
velocity), it proves the exponential asymptotic decay of a solution of the linearized 
Vlasov equation in a bounded spatial domain, as predicted by Landau. However, it is 
known that such a linear solution may decay at a much slower rate, and even not decay 
at all, when the analyticity hypothesis for the perturbation (see [U [5] for references 
in the physics literature) or the bounded spatial domain hypothesis [6] is not satisfied. 
Mouhot and Villani proved recently the asymptotic exponential decay of a perturbation 
evolving according to the full non linear equation in a bounded domain, using analytic 
norms (which implies the analyticity of the stationary state and its perturbation) (TJE]- 
However, it is shown in [9] (see also [10] for a previous non rigorous treatment) that 
such a non linear damping may fail if weaker norms are used. 

In many cases, one would like to ask the same question in the case of non 
homogeneous stationary solutions. In particular, this is the case in the astrophysical 
context, where the concept of Landau damping is often used [TTJ. [2] contains a 
discussion of the linearized Vlasov equation around a non homogeneous stationary 
solution. However, the analysis of this situation faces some technical difficulties. At a 
formal level, these difficulties may be partially overcome using the "matrix" formulation 
introduced in the works of Kalnajs [12] and Polyachenko and Schukhman [13]. Since 
then, many papers have computed linear instability rates in an astrophysical context 
using this method (see [HI [15], to mention just a few). Recently, some new methods to 
investigate the stability and compute such growing rates for unstable non homogeneous 
stationary states have been introduced and tested on toy models [HI [TTJ [HI [H] . Beyond 
these unstable states, stable oscillating modes in the context of ID self gravitating 
models have also been investigated in some details [201 E]. However, there seems to 
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be very few analytical studies of what would be the strict analog of Landau damping: 
the decay of a perturbation close to a stable non homogeneous solution of the Vlasov 
equation. The reason is probably that such a study is technically more difficult, since 
it requires the further step of an analytical continuation. Weinberg has performed this 
step using a numerical approximation and computed the analog of a Landau damping 
rate in some cases (22] (see also [23]). However, at variance with the homogeneous case, 
we shall see that this Landau damping rate never controls the asymptotic decay of the 
perturbation. 

It is well known that the 2D Euler equation, as well as other related conservative 
2D fluid equations, share a lot of similarities with the Vlasov equation. There have 
been many linear investigations of perturbations around stationary base flows such as 
a vortex or a shear flow, starting with Rayleigh [24]. Transposing these studies in 
the context of the Vlasov equation, linearized around a non homogeneous stationary 
state, one would expect the following generic picture: the dispersion relation has branch 
point singularities on the real axis, which make the analytical continuation procedure 
used in Landau's work trickier; if the stationary state is stable, the asymptotic decay 
of the perturbation is controlled by these branch point singularities, and is in general 
algebraic. An exponential decay of the perturbation "a la Landau" may be visible, 
but it is restricted to an intermediate time window, before the algebraic decay kicks 
in. In particular, it seems reasonable to expect that the exponential damping studied 
by Weinberg in the gravitational case [22] is actually a transient effect (this statement 
does not preclude of course its potential physical relevance). We remark that such an 
exponential decay followed by an algebraic decay has also been studied in a dissipative 
system of coupled oscillators [25] ; the decay mechanism in this case is similar to the one 
in Vlasov systems. 

Smereka [26] investigated the asymptotic behavior of a linearized ID Vlasov 
equation around non homogeneous stationary states, and reached the conclusions 
outlined above; his analysis is however limited to a particular class of interactions, 
and does not contain direct numerical simulations. Moreover, his conclusion is affected 
by a non generic choice of the perturbations, as we will discuss later on. Besides [26J, 
we are not aware of other studies tackling this problem. 

To be more specific, our goals in this paper are 

(i) Study as generally as possible the linear asymptotic decay of perturbations around 
non homogeneous stationary states of a ID Vlasov equation. 

(ii) Perform explicitly the computations in the simple case of the Hamiltonian mean- 
field (HMF) model. 

(iii) Compare the analytical results with detailed numerical computations; this is made 
possible by the use of a simple toy model such as HMF. 

This paper is organized as follows. All theoretical results are described in SecEJ 
From the linearized Vlasov equation, we derive two equations for the Fourier-Laplace 
components of perturbation and of potential in Sec J2.ll In Sec l2.2l these two equations 
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are formally solved in the Laplace space with the help of biorthogonal functions; an 
example of such functions is given in Sec J2.3l for spatially periodic systems. The 
asymptotic dynamics of a perturbation is determined by its singularities in the Laplace 
space; we classify the singularities in Sec J2.4l Focusing on one type of singularities, 
which is in many cases the relevant one, we show that the perturbation asymptotically 
decays algebraically with the exponent —2 in Sec J2.5l These general results are applied 
to the HMF model in Sec J2.6| in this case, we show that the two components of the 
magnetization vector decay with exponents —3 and —2 respectively. We note that the 
exponent —3 comes from a special cancellation due to a symmetry of the HMF model. 
The theoretical results for the HMF model are numerically examined in Secj3j In order 
to test the two exponents separately, we introduce two types of perturbation in Sec J3.ll 
and in Sec 13. 21 respectively. The last section |4] is devoted to conclusion. 

2. Theory 

We consider the Vlasov equation in one dimension for the one-particle distribution 
function f(x,p, t), 

d t f + d p Hd x f - d x Hd p f = 0, (1) 

where x G D C R is the position variable, p G R the conjugate momentum variable, 
and H is the one-particle Hamiltonian defined by 

H[f](x,p,t) = ^ + <j>\J\{x,t) + cxt (x). (2) 

The potential (f)[f](x,t) is defined from the two-body interaction potential v and the 
distribution / as 

<!>[f}{x,t)= v(x -y)f(y,p,t)dpdy, (3) 
J d Jm 

and (pex.t(x) derives from an external force F ext (x) as 

^extO) = -d x (j) cx t(x). (4) 

Although it would be interesting to consider a time dependent and/or non-potential 
external force, we focus in this article on a static and potential force such as (j4"l). The 
domain D is typically R or [0, 2n}: 

(i) Case 1, D = R: The model is defined on the whole real line. A typical example is 
a ID self-gravitating system, a very much studied caricature of the more realistic 
3D self-gravitating systems. 

(ii) Case 2, D = [0, 2ir]: The system has periodic boundary conditions. Such boundary 
conditions are sometimes used in plasma physics; it is also the setting of the HMF 
model, a paradigmatic toy model for long range interacting systems (see Sec. I2.6j) . 

We give in the following subsections a general analysis of the linearized Vlasov equation, 
which applies to the both cases. 
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2.1. Two equations to be solved 

Let us call fo a stationary solution of the Vlasov equation (DO). The one-particle 
Hamiltonian H[f Q }(x,p) is integrable, and its trajectories are level curves of H[f }. If 
D = R, one could imagine trajectories that are unbounded and not periodic. However, 
the stationarity of fo imposes that fo(x, p) is constant along the trajectories; this imposes 
that such unbounded non periodic trajectories either do not exist, or are not populated 
by the density fo. As a consequence, we can always define the angle 9 and the action 
J from the original coordinate (x,p). Strictly speaking, this change of variables is not 
always one-to-one: whenever there exists a separatrix in the one-particle Hamiltonian 
H[fo](x,p), there are distinct trajectories with the same value of the action. Thus, 
although we will formally use this change of variables, a careful treatment may be 
needed for specific potentials (see for instance [28]). 

We add a perturbation /i to the stationary solution, and start from an initial 
condition fo + fi(t = 0). The Hamiltonian H is linear with respect to /, and we have 
H\f + f 1 ]=H[fo] + (f>\fi], where 



= 0iOM) = / / v(x - y)ft(y,p,t)dpdy. (5) 




' D 



The stationary solution fo must be constant along trajectories of the Hamiltonian H[fo\; 
a sufficient condition for this is to take fo a function of the action alone. It is not a 
necessary condition if there is a separatrix, since two disjoint trajectories correspond to 
the same action. Thus, with a slight loss of generality, we will assume in the following 
that fo may be written as fo{J)- The unperturbed part of one-particle Hamiltonian 
is also a function of the action alone, and written as H[fo](J). Using the angle-action 
variables (6, J), we have the linearized Vlasov equation 

ft/i + n(Wi-/ (Wi = o, (6) 

where we have defined the frequency Q(J) = dH[fo]/dJ. Notice that the external 
potential does not enter in this linear equation; it appears implicitly of course through 
the definition of the angle-action variables. 

To analyze the linearized Vlasov equation ([6]), we introduce the Fourier-Laplace 
transform u(m, J, uf) of a function u(9, J, t) as 

d9 e~ imd / dt e iujt u(6,J,t) (7) 

■IT JO 

where m is an integer and Im(o;) large enough to ensure convergence. The inverse 
transform is then 

1 °° f 

u(6, J, t) = I duj ^( m ' J ' u)e~ luJt e ime (8) 



r 



where T is a Bromwich contour running from — oo + io to +oo + io~, and the real value a 
is larger than the imaginary part of any singularity of u(m, J, u) in the complex w-plane. 
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Performing a Fourier transform with respect to 9 and a Laplace transform with 
respect to time on ©, we obtain, after simple algebraic manipulations 

fi(m, J, oj) = A(m, J, u)<fii(m, J, u) + B(m, J, uj), (9) 

where 

A(m,J,«,)= ™f° T [ J) , (10) 

nt t \ dim, J) , s 

B(m,J,u)= , 11 

m\l{J) — uj 

and ig(m,J) is the Fourier transform of the initial perturbation fi(9,J,t = 0) with 
respect to 9. We assume 

J J f 1 (9,J,t = 0)d9dJ = (12) 

and hence g(0, J) = 0. 

The two equations (jSJ) and relate fx to <f>\. The strategy is now to combine these 
two equations to compute f\ and <f>\. One sees however that fx is easily obtained in 
angle-action variables whereas <f>x is more easily expressed in the original (x,p) variables. 
To overcome the difficulty of the two natural coordinate basis (x,p) and (9, J), we follow 
the standard procedure and introduce two families of biorthogonal functions [271 E21 EES] . 

2.2. Biorthogonal functions 

We introduce the linear mapping L v by 

L v : di-> u (13) 
where u = L v ■ d is defined by 

u(x) = v(x- y)d(y)dy. (14) 

J D 

We assume that there exist two index sets I and I' which satisfy I' C / C Z, and 
two families {ci,-(x)}j e / and {«&(£)}-£<=// which satisfy the following conditions: 

(i) {dj(x)}j£i is linearly independent and any density function p(x) may be expanded 

as 

P( x ) = ^2 a j d j( x ), ( 15 ) 
i ei 

(ii) {uk(x)}kei> is linearly independent and spans Im(L„); any function (?(x) G Im(L„) 
may be expanded as 

9(x) =^2b j u j (x), (16) 

(iii) the two families are orthogonal to each other: 

(dj,u k ) = / dj(x)u k (x)dx = \kd~jk, (j € I,k e I') (17) 
J D 

with Afc 0, and where 5jk is the Kronecker S. 
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(iv) For all k G I', dk and Uk satisfy the property Uk = L v ■ dk- For all k G I \ F, 
L v ■ d k = 0. 

Let us start from the definition of the potential <fri(x, t) ([5]). This definition may be 
rewritten as 



0i(x,t)= / v(x - y)pi(y,t)dy (18) 

J D 

by using the perturbation density 

Pi{x,t) = / f 1 (x,p,t)dp. (19) 

From the assumption (i), the perturbation density is expanded in the form 

pi (z, t) = ^ a i (*) d i ( x ) • ( 20 ) 
j el 

Substituting (1201) into (TTBl and using (iv), we obtain an expansion for X in the form 

0i(x, t) = ^ a k {t)u k {x). (21) 
fee/' 

The Fourier-Laplace transform of (121]) is expressed by 

0i(m, J, w) = ^ a fe (cj)c fcm ( J), (22) 

where dk{oj) is the Laplace transform of ak{t), 

POO 

~a k (u)= / a k (t)e iujt dt, (23) 
and Ckm(J) is the Fourier transform of Uk(x), 

CkM) = r u k (x)e- me d9. (24) 



Substituting ( 1221) into (Ej), we obtain 

/i(m, J,w) = A(m, J, u) ^a fc (u;)c fcm (J) + B(m, J, a;). (25) 

fee/' 

The inverse Fourier transform of ( l25l) gives fi(8,J,u), the Laplace transform of 

A(0, J, a,) = ±- V A(m, J,u)e imd (26) 

/7T g ■* 



-E 

2tt ^ 



A(m, J, cu) ^ a k (uj)ckm(J) + JB(m, J, oj 

fee/' 



im8 



(27) 



We observe that J, cu), and hence fi(9,J,t), are determined by the {5fc(u;)}fe e j/, 
so that the subset {%(w)}je/\/' is not necessary. We therefore seek a solution for the 
subfamily {ak{w)}ker instead of the whole family {aj(cu)}j e i. For this purpose, we 
multiply fl26|) by n^(x) (/ G i 7 ) and integrate it over 9 and J. Using the biorthogonality 
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relation (j!7p and noting that the change of variable (x,p) h-> (0, J) is symplectic 
dx A dp = d9 A dJ, the left-hand-side of fl26l) becomes 





fi(8, J,u)ui(x)d9dJ = J J fi(x,p,u)ui(x)dxdp (28) 

Pi(x, oj)ui(x)dx = ai(oj)\i. (29) 



D 



To derive the last equality, we have used the fact that the Laplace transform of pi(x,t) 
dm is 



Pi(x,u) = a,j(u)dj(x). (30) 

On the other hand, the right-hand-side of (J26J) , submitted to the same operations, 
becomes 

J J —^^A(m,J,cu)^2a k (cu)c km (J)+B(m,J,cu) ui(x)e imd d9dJ 

J J mGZ fee/' 

= 2~ Z Z / ^^ m ' J ' U ) C km(J)ci m (J)dJ 

mezkei' 

I B{m,J,uj)c lm {J)dJ. (31) 
Remembering the definitions of A(m, J, oo) and B(m, J,u), we introduce the functions 
Flk ^ = J" Z) / m fl/n_ c "m(J)ci m (J)dJ, l,kel' (32) 

27r / miz J — a; 

where contributions from m = vanish not only for Fik but also for Gi, thanks to 
assumption ffl2|) . We further define the (jji 7 ) x (jji 7 ) matrices A and F(u>) = (i 7 /fe(w))/,fee/' ) 
where A is diagonal with elements {Xk}ker, and the (jj/')-dimensional vectors G(u) = 
(Gi(u))i e p and a(oj) = (ai-(uj))kei'- Using the above matrices and vectors, the equation 
for a(oj) reads in matrix form: 

[A-F(u)]a(u) = G(u). (34) 

The equation ( 1341) is formally solved as 

a{u) = [A-F(u)]- 1 G(u)), (35) 

and the temporal evolution of the {a,k(t)}kei', and of fi, is obtained from the inverse 
Laplace transform of d(u). The inverse matrix [A — F(u)}~ 1 does not always exist since 
the determinant det(A — F(u)) is not always non-zero. This determinant is sometimes 
called the dispersion function, and its roots are poles of d(u). We will discuss the 
singularities of a(uj) in Sec J2.4l after giving an example of the two families {dj}j £ i and 
{linker in the domain D = [0, 2tt]. 



Algebraic damping in the one- dimensional Vlasov equation 9 
2.3. Example: D = [0,2tt] 

If we consider the domain D = [0, 2tt] with periodic boundary condition, the interaction 
potential v(x) must be also 27r-periodic and is expanded in Fourier series as 

= ;r Z>™ eira ". ( 36 ) 

where the coefficients determined by 

v(x)e~ imx dx. (37) 



We can choose the two families {dj(x)}j^i and {uk(x)}kei> as 

d J (x) = ^ x , jel = Z, (38) 

and 

p2tt j 

u k (x) = v(x- y)d k (y)dy = —v k e tkx , k e I' (39) 
Jo 271 

where the index set I' is defined by 

/' = {k e z | v k ^ o}. (40) 

From the definition of Uk{x) and I' the assumption (iv) is satisfied. Expansions on the 
two families {dj(x)} j e z and {uk(x)}k£p are essentially Fourier expansions, and hence 
the assumptions (i)-(iii) are also satisfied. The factors are = Vk/27r for k G /'. 
Generically, all are non zero, and I = V = 7h. However, for the HMF model, which we 
will introduce in section [2161 v{x) = —cosx and lm(L v ) is 2-dimensional. Accordingly 
the index set I' is I' = {1, — 1} and the matrix A — F(u) in (|34|) is a 2 x 2 matrix. 



2.4- Singularities of a{uj) 

From the knowledge of the functions a(u), one may easily compute the time evolution 
of the potential and density perturbations, through an inverse Laplace transform. The 
asymptotic in time behavior of this inverse Laplace transform will be determined by the 
singularities of the functions dk{oj) (k G I') in the complex plane. We turn now to the 
study of these singularities. 

The matrix coefficients Fik(cu) and the vector coefficients Gi(u) are defined through 
integrals over the real variable J, see (132]) and fl33|) . These integrals are naturally defined 
in the whole half plane lm(u) > 0. Note however, that expressions ( 1321) and ( |33j) are 
in general not properly defined for u G R, as in this case mVt(J) — u may vanish. For 
lm(u) < 0, we will actually have to consider rather the analytical continuations of the 
expressions fl32l) and fl33|) . We may have two kinds of singularities, described in the 
following. 

The first kind of singularities is poles, coming from roots of the dispersion function 
det(A — F(cj)). If such a root exists in the half plane Im(cj) > 0, it corresponds to 
an eigenvalue of the linearized Vlasov operator, and it yields an exponential growth of 
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the perturbation. In the upper half plane Im(oj) > 0, the only possible singularities for 
a (a;) are the poles. We assume in the following that the reference stationary state /o 
is linearly stable, so that the determinant of A — F(u) does not have any roots with 
lm(u) > 0. If det(A — F(oj)) has a root on the real axis lm(u) = 0, it corresponds to a 
purely oscillating mode. This is compatible with a linearly stable fo, but corresponds to 
a non decaying perturbation; we also assume in the following that this does not happen. 
The analytical continuation of det(A — F(u)) in the lower half plane Im(a;) < may also 
have roots. They correspond to "Landau poles" for a(u), and give rise to an exponential 
damping of the perturbation. This damping behavior is known as Landau damping, and 
has been studied in the gravitational case in [22} [23], and in the HMF case in [283 . As 
anticipated in the introduction, we will see that this exponential damping, if it exists, 
is subdominant in the large time region. 

The second kind of singularities comes from the integral in ( 132]) and ( 1331) . and 
appears on the real axis of Im(co>) = 0. We have to study the singularities of functions 



properly defined for Im(2;) > 0, where ip and Q are real functions, and [a, b] C RU{+oo}. 
J has to be thought of as the action coordinate, and fl as the associated frequency. We 
assume that ip is analytic. We may set m ^ since the contributions from m = in 
(13"2"|) and f l3"3"|) vanish. We now show that tp(z) is regular for z G R except for special 
points, and will classify the special points into three types. Notice that both functions 
Fki and Gk fit in this framework. 

If the equation mO(J) = z has no solutions in [a, b] for any z in a neighborhood 
of real zq, then if is analytic in a neighborhood of zq. Assume now that the equation 
mQ(J) = z has one or several branches of solutions J*(z) G [a, b] in a neighborhood 
of z , where all J* (z) G [a, b] are regular as functions of z. In this case, <p can be 
analytically continued from the open half plane Im(2;) > to the neighborhood of 
z by taking into account the possible residue contributions of the roots J*(z), in 
a straightforward generalization of the "Landau prescription". Thus, generically, no 
singularity of ip appears at z = z . 

The singularities of <p are hence associated with special points zq, such that the 
branches of solutions J*(z) G [a, b] of the equation mQ(J) = z undergo a bifurcation or 
are singular. This may happen in the following three types, illustrated on Fig. [Q 

(i) zq = mfi(fl) or zq = mQ(b), when one or both are finite. This is a common 
situation, generically encountered in ID self-gravitating systems [20] as well as in 
the HMF model [28]. For instance, this mechanism creates a singularity around 
the frequency muo = mQ(J = 0), where J = corresponds to the minimum of 
the effective potential 0[/o] + </>ext> see Fig- CD We will show in Sec J2.5l that this 
singularity is logarithmic. 

(ii) Zq corresponds to a J such that Q(J) is singular: this may be an action J s 
corresponding to a separatrix; an illustration is given on Fig. [fl for J = J s , u = 0. 




(41) 




Figure 1. An example of function 0(J), with the special points at the origin of 
singularities for the associated function ip defined as in Eq. (|4Tj) . 



This is not a common situation for the ID self-gravitating models we have in mind. 
This generically happens however for periodic systems, where the trajectories in the 
one-particle Hamiltonian H[f ] may be oscillating or librating, and the two regions 
are delimited by a separatrix. We will see such a situation in the HMF case. 

(hi) Zq corresponds to a J such that VL (J) = 0: this corresponds generically to a local 
maximum or minimum of the frequency (see frequencies u\ and U2 on Fig. [TJ). 
We are not aware of any model studied in the literature where this phenomenon 
happens. This is certainly an interesting case to study, especially in view of 
the results obtained in an analogous situation for the diocotron instability of a 
magnetically confined electron column [29] . and recently for the 2D Euler equation 
[30] . We will not be concerned with this type of singularities in the following. 

In the next subsection we estimate the asymptotic relaxation of a>k(t) by considering 
contributions from singularities of the first type, since they should be in many cases the 
relevant singularities. We will confirm whether the estimation is valid by performing 
direct iV-body simulations in Secj3j 

2. 5. Contribution from singularities of the first type 

We concentrate now on the singularities of the first type at z = mfl(a) = mu a . Notice 
that if b = +00, and = fi(+oo) is finite, z = mcu^ is a singular point; indeed, the 
number of solutions of the equation mQ(J) = z changes from z > mw M to z < rriUoo. 
We expect this singularity to be irrelevant if fo(J) decreases rapidly enough as J 00, 
and we neglect it in the following. 

To analyze the function <f(z) fj4"T|) around z = mu a , we expand Q(J) and ip(J) in 
power series around J = a as 

Q(J) = u a + c(J - a) + O ((J - a) 2 ) (42) 

and 



1){J) = d(J-aY + o((J-a) u ) 



(43) 
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for some v G N. The exponent v will be determined for a given m later. Substituting 
these assumptions into (p and changing J — a to J, this yields the following singular part 
for ip: 

rb—a jv 

ip(z) = C - — -dJ + (regular part) (44) 
Jo -J ~ C 

with the constant C = d/mc and ( = (z — mu a )/mc. Using the equality 

TV JV-1 

r- 1 + C^- (45) 



j-C v-c 

recursively, the singular part of <p(z), which comes from the lower bound of the integral 
( 1441) . behaves clS 8b logarithm times a power around £ = 0: 

(f(z) = cste x (z — mu a y \\i{z — mu a ) + (regular part). (46) 

Thus, at u = mu a , the matrix elements Fki(u) and the vector elements Gk{oo) have 
singularities of the type (u — mu a ) v 'ln(w — muj a ). From (|35|) . dkiui) is expressed as a 
sum, product and ratio of Fki(uj) and Gk{oj) functions. Thus, the leading singularity of 
akijjj) at oj — mojQ is also of the type (u — mu) a ) v 'ln(a; — muj a ) for some z/. afc(t) is the 
inverse Laplace transform of 0^(0;). Since a& has only logarithmic singularities on the 
real axis, we may deform the contour V of the inverse Laplace transform down to the 
real axis. The inverse Laplace transform then becomes an inverse Fourier transform. 

The asymptotic decay of a^(t) is then determined by the strongest singularity of 
ak{uj) on the real axis (see [31] p. 52). A singularity such as PS"]) yields an asymptotic 
decay as (see [31] p. 42) 

g— imujat 

a k (t) ~ cste-^^. (47) 

See Appendix A| for a heuristic explanation on how to obtain estimates such as fT4T|) . 

To obtain v for a given m, we need to introduce some assumptions about the system 
we consider. We assume that the potential created by the stationary state (interaction 
potential + external potential (j)[f }(x) + e xt(^)) has a single minimum at J = a = 0, 
and is quadratic with respect to x around its minimum at leading order, except for 
an irrelevant constant term. This is the case in many situations of interests, such as 
stationary states for a ID self-gravitating system. We also assume that /o(<^) is analytic, 
and decays fast enough at infinity (for instance exponentially). This assumption excludes 
for instance truncated /o ; or compactly supported stationary states. As an example, 
the thermal equilibria of a ID self-gravitating system or of the HMF model satisfy all 
assumptions. Finally, we assume that the perturbation is also analytic. 

Under the above assumptions, we now estimate the exponent v. From (132]) and ( 1331) . 
we see that the function ip(J) reads 

^{J) = fo(J)ckm(J)cim(J) 
for Fki{u) functions and 

ip(J) = g(m, J)c lm {J) 
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for G(oj) functions; we now expand these functions with respect to J. Let us start from 
the function Cfc m (J). In the limit of J — > 0, motion is harmonic with the frequency 
uq, and the position x is written in polar coordinates using angle- action variables as 
x oc J 1 / 2 sin 9. Expanding Uk(x) with respect to x, and substituting the above expression 
of x into the expansion, the function Cfc m ( J) reads 

"2tt 00 n ,(n)i 



^° n=0 



where iz^ is the n-th derivative of Uk- The first non- vanishing term corresponds to 
n — \m\, and hence the leading order for Ck m {J) in a small J expansion is 

c km {J) ~ J |m|/2 - (48) 
Ckm(J) is the Fourier transform of Uk{x); in a similar way, the leading order in a small 
J expansion for g(m, J), the Fourier transform of the initial perturbation, is 

g{m,J)~j\ m \l\ (49) 

since the perturbation is regular. The function f' Q {J) is regular and hence the leading 
order is constant. Consequently, the function ip(J) is, at leading order 

ij){J) ~ J |m| (50) 

both for F and G functions. Hence, for each term in the infinite series defining the 
coefficient Fkiius) ( |32l) and Gi{uj) ( l33l) . we have v = \m\. 

Going back to (j3"21) , and making use of section 12.41 we see that the strongest 
singularities for Fj,i{lo) and Gi(u) come from the m = ±1 terms in the sum over m, so 
that the exponent is v — 1. We conclude using (T4T1) that the functions afc(t), under the 
hypothesis of this section, decay as e~ lulot /t 2 . 

This result has to be compared with [26], which finds a decay exponent 3/2. 
Since we have performed the same kind of analysis as Smereka does in [26], this 
discrepancy is surprising even if the class of Hamiltonians studied is different. It 
may be traced back to the fact that this author uses the following hypothesis for the 
perturbation: \im.j^ g(m, J) ^ even for m/fl; this would correspond to a singular 
perturbation, since the initial perturbation fi is not well-defined in the limit J — > 0, 
because limj_;. e me g(m, J) ^ limj_> e ime> g(m, J) as soon as e im< y e - e ') l. This singular 
perturbation implies that g(m, J) ~ J° in the limit J — > instead of ( )49l) . and hence 
ijj(J) = g(m, J)ci m (J) ~ jH/ 2 . Accordingly, the strongest singularities for G[(cu) ( 1331) . 
coming from the m = ±1 terms, corresponds to v — 1/2. Using considerations similar 
to the ones described in Appendix A Smereka showed that this gives a decay exponent 
3/2. Considering in Smereka's setting a regular perturbation, as is more natural (and 
as actually assumed in Eq. (21) of [26]), would yield the —2 exponent also for the class 
of Hamiltonians studied in [26J . 

We may also compare this result to the asymptotic decay of perturbations around a 
stationary 2D shear flow or a vortex. In these cases, the longitudinal (resp. transverse) 
velocity perturbation asymptotically decays as 1/t (resp. as 1/t 2 ), without temporal 
oscillations. 
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2. 6. Example of the HMF model 

In this section, we analyze in more details a specific example, the HMF model whose 
Hamiltonian is 

N 2 -y N N 

H(x, P ) = E| + ^EE[ 1 - cos ^ - **)]■ ( 51 ) 

j=l j=l k=l 

The canonical equation of motion of the HMF model is described through the 
magnetization M = (M x ,M y ) defined by 

1 N 1 N 

Mx = ]yYl cos x 3> M v = nY1 sin X J' ( 52 ) 

3=1 3=1 

and hence the computational cost is O(N) for each time step, although the number of 
interactions between the N particles is 0(N 2 ). This advantage allows precise numerical 
tests of the predictions. The associated Vlasov equation reads: 

dl +p d_l_dj[f]d_l = ^ 

dt dx dx dp 



with 



and 



(j>[f]{x,t) = -M x [f] cosx - M y [f] smx (54) 



/OO /*27T 
/ cosx t)dxdp, (55) 
-oo Jo 

/oo p'2ti 
\ smxf(x,p,t)dxdp. (56) 
oo JO 

Note that this is a pendulum potential, so that the dynamics admits a separatrix. The 
action-angle variables are explicitly written in terms of elliptic integrals [28]. Without 
loss of generality, we consider a stationary solution with M y = 0, and write Mx (t) and 
My 1 " 1 (t) the magnetization perturbations. 

For convenience, we choose real functions for the family {uk}kei r - u c (x) = cosx and 
u s (x) = smx and /' = {c, s}. The coefficients a c (t) and a s (t) of the potential <f>i(x,t), 

4>i(x, t) = —a c (t) cosx — a s (t) smx, 

correspond to M^\t) and My (t) respectively. 

The expansions of u c and u s in the Fourier series of the angle variable 9 define 
the coefficients {c cm (J)} mg z and {c sm (J)} me z as in (124]) ; to simplify the notations, we 
rename these coefficients c m (J) and s m (J) respectively. 

This choice for the family {uk} makes the matrix A — F(u) diagonal [28], with 

Fcc(u) = f £ / on[ J) MJ)\ 2 dJ, (57) 
2tt ^— ' / mill J) — oo 

Fss(uj) = f Yl I nn\ J) \sm(J)\ 2 dJ, (58) 
2tt z — ' / mil (J) — oo 

F cs (oo) = F sc (oo) = 0. (59) 
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We also have the following expressions for G: 

G -M = f (60) 

m£Z J y ' 

i y r 3M ( J)dJ (61) 



^4 priori, according to the general discussion in the previous section, the leading 
singularity of F cc and F ss is located at uj = uq = y/Mo, and has an index v — 1, since 
Wo here plays the role of uj a in the previous section 12.51 However, due to the symmetries 
of the system, a further cancellation occurs: the function c m (J) identically vanishes for 
all m odd and all J < J s , with J s the action at the separatrix. Thus, the strongest 
singularity for F cc and G c actually comes from the m = 2 term in ( ISTjl and ( )60l) . and 
is located at a; = 2wo; its index is v = 2. This has an interesting consequence on the 
asymptotic behavior of M^\t)\ since it is now governed by a singularity with index 
v = 2, we expect 

e -2ioj t 

^fort^oo (62) 

at variance with My l \t), which is still governed by a v = 1 singularity 

p—iuot 

M^\t) — fort^oo. (63) 

This feature makes the HMF model particularly suitable for a numerical test of the 
theory developed in this section, as we should be able to probe two different asymptotic 
behaviors for M { }\t) and M^\t). 

3. Numerical simulations 

In this section we test numerically the linear predictions of the previous section, on 
the example of the HMF model, by solving the whole (non linear) Vlasov equation. 
Solving the Vlasov equation over long times may be a very heavy numerical task, or 
even impossible with current computers. In this case, several features help: the model is 
one dimensional, it is particularly simple, and we only need to solve the Vlasov equation 
close to a stationary state. 

A natural strategy could be to solve directly the Hamiltonian iV-body dynamics, 
with N large enough; we know that this provides an approximation to the continuous 
Vlasov evolution. We have found that the fmite-iV fluctuations were too big to 
allow a test of the asymptotic in time regime. Another strategy would be to use a 
standard Vlasov solver; for instance, a semi-Lagrangian method has already been used 
for HMF [32]. This resulted in very heavy computations. Finally, we have chosen to 
introduce an algorithm relying on a Hamiltonian simulation of appropriately weighted 
particles [33] , it provides a very convenient tool to test the theoretical predictions. The 
algorithm is described and discussed in Appendix B| We have tested the weighted 



particles algorithm against (i) a semi-lagrangian code (ii) a simple unweighted iV-body 
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code (results not reported). The temporal evolution of the magnetization from the three 
codes are in good agreement up to a certain time. The unweighted iV-body simulation 
becomes dominated by finite size fluctuations much earlier than the weighted particles' 
one. 

All simulations discussed in the following were performed using the weighted 
particles algorithm, close to a thermal equilibrium stationary state, parametrized by 
the temperature T = 1//3: 

f (x,p) =Ue^ p2/2 ~ MoCOSX ) (64) 

where M is the normalization and the magnetization M (T) is solution of a consistency 
equation [M] . The thermodynamical equilibrium state of the HMF is non homogeneous 
(that is Mq 7^ 0) as soon as T < 0.5. In the following we only use T = 0.1; larger 
temperatures resulted in increased fluctuations and made it more difficult to reach 
the asymptotic in time regime. The magnetization is M = 0.946 and the harmonic 
frequency is ujq = y/Mo = 0.972 for T = 0.1. 



3.1. Cosine perturbation 

We consider first a cosine perturbation of the thermal equilibrium: 

fo(x,p) + h{x,p) = A4 e -Kp 2 /2-Mocosx) (1 + acogx) ^ (65) 

where \a\ is small enough. This perturbation is compatible with the (x,p) — > (—x, —p) 
symmetry of the canonical equation of motion, and (t) is then identically equal to 
zero. We may restrict the initial points such that Pi(0) > 0, and obtain the temporal 
evolutions of particles which are initially in the lower half of /z space without direct 
computations, using this symmetry. 

To estimate the perturbed magnetization M^~\t), we subtract from M x (t) its long 
time average. A typical temporal evolution of M^\t) is shown in Figj2j Finite size 
effects, which are visible on the curve for N = 10 7 points, may prevent the study of the 
asymptotic behavior, so that it is usually necessary to use a very large number of points 
(see the curve for N = 10 8 points). We fit the envelop of the decaying M^\t) curve 
by a power-law, using the least square method (see figure Q). We find an exponent 
—2.95, which is in very good agreement with the prediction (162"]) . 



3.2. Sine perturbation 

We consider now a sine perturbation of the thermal equilibrium: 

f (x,p) + h{x,p) = A/^-K^-Mocos*) ^ + aginx) _ (g6) 

In this case, the (x,p) — > (—x, —p) symmetry is broken, and we have to compute the 
evolution of the whole /i space. A slow rotating motion of the magnetization appears, 
which makes it difficult to define an asymptotic average value of Mi l \t) and M^\t). 
Rather, we used a running average to eliminate the rotation effect (see Appendix Cp . 
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semi- Vlasov code (n=w 7 ) 
semi- Vlasov code (jV=lo rt ) 



Figure 2. Temporal evolution of M^\t) for an initial condition as in (|65|) . We used 
the weighted particles code with T = 0.1, a = 0.1, N = 10 7 (black curve) and N = 10 s 
(gray curve). The initial weighted points are equally distributed in ] — tt,it} x [—3,3], 
and we used the ix,p) — > (— x, —p) symmetry in order to reduce the computations. 




Figure 3. (color online) Temporal evolution of the M^~\t) envelop for an initial 
condition perturbed by a cosine, as in (j65[) . Red points: numerical simulation using 
the weighted particles code with T = 0.1, a = 0.1 and N = 10 9 . The initial weighted 
points are equally distributed in ] — tt, it] x [—3, 3]. Black dashed line: power law fitting 
from t = 600 to 6000 (oc i" 2 - 95 ). 



Finally, we obtain curves similar to the case of a cosine perturbation to confirm the 
prediction ( |62|) and ( )63l) . 

We then fit the envelops of M^\t) and M^\t) with power laws. The exponents 
are —2.94 and —1.77, to be compared with the predicted —3 and —2, see ( 1621) and ( 1631) . 
The relative error for M^\t) exponent is close to 10%. Different explanations are 
possible: the fit range does not completely lie in the asymptotic regime and/or there are 
numerical errors. ( 162|) and (|63l) predict an oscillating decay of the perturbations (t) 
and My(t), with frequency respectively u and 2tu - On Fig. [51 we plot the power 
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t 

Figure 4. (color online) Temporal evolutions of the M^\t) (upper panel) and My 1 ^) 
(lower panel) envelops for an initial condition perturbed by a sine, as in (|66j) . Red 
points: numerical simulation using the weighted particles code with T = 0.1, a = 0.1 
and N = 10 9 . The initial weighted points are equally distributed in ] — tt, it] x [—3, 3]. 
Black dashed line: power law fit from t = 600 to 6000. We obtain an exponent equal 
to -2.94 for Mi 1] (t) and -1.77 for M^\t). 

spectra of M^\t) and My(t). We observe that indeed in both cases a single frequency- 
is selected in the long time regime, with numerical values respectively 1.944 and 0.973. 
This is in almost perfect agreement with the theoretical prediction 2uq = 1.945 and 
u = 0.972. 

Figure [5] also allows to observe the cross-over between two different dynamics 
explored by the system: the short time evolution is driven by the Landau pole 
contribution, and the asymptotic behavior by (|62|) and ( )63|) . Indeed, as shown in [28] . 
it is possible to compute the dominant Landau pole for the parameters of Fig. |5j one 
finds a frequency Re(cjL) — 1-8. This is in good agreement with the short time power 
spectrum of M x (t), which is maximal around u = 1.85. We therefore conclude that 
Landau damping occurs in the short time regime, before the algebraic decay dominates. 
The peak position at short time in the M y 's spectrum may be a signature of the growing 
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ojq peak in a long time regime, since it is close to ojq. It might also be related to a Landau 
pole associated to M y : since the matrix A — F(oj) is diagonal, Landau poles for M x and 
My may be different. Finally notice that the power spectrum divergence close to u = 
is only due to the rotation of the magnetization. 



DFT of M,(() for < f < 200 o- o 
DFT of M x {t) for 600 < t < 6000 •— i - 
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Figure 5. (color online) Power spectra of the M^\t) (upper panel) and My X \t) 
(lower panel) for an initial condition perturbed by a sine, as in (|66[) . The numerical 
simulation was done using the weighted particles code with T = 0.1, a = 0.1 and 
N = 10 9 . The initial weighted points are equally distributed in ] — it, ir] x [— 3, 3] . Each 
inset includes a magnification around the maxima and the short time power spectrum. 
The maxima of the long time power spectra of (t) and M^p (t) are respectively 
1.944 and 0.973. 



4. Conclusion 

We have investigated the asymptotic dynamics of perturbations around stable non 
homogeneous backgrounds in spatially one-dimensional Vlasov equations. The 
dispersion relation of the linearized Vlasov equation has poles in the lower half of 
the complex plane and logarithmic branch points on the real axis. The poles yield 
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exponentially decaying contributions: this is a form of Landau damping. The branch 
points yield algebraically decaying contributions. If the dominant branch point is related 
to a minimum of the potential created by the stationary state, the leading order of 
the potential is quadratic with respect to position around the minimum, and unless 
some special cancellation occurs, the perturbation potential behaves asymptotically as 
e™ 0t /t 2 , where Uq is the harmonic frequency of the potential well. We expect this 
situation to happen in many cases of interest. 

We have tested the theory on the HMF model by performing iV-body simulations 
which correspond to the full Vlasov equation in the large iV limit; these simulations 
used the weighted particles code. The exponent and frequency of the decay have been 
confirmed by these simulations, including in one case where a special symmetry imposes 
a decay as e i2uJot /t 3 . 

This summarizes in the following scenario for decaying perturbations around stable 
non homogeneous background in one-dimensional systems: 

i) The perturbation potential first roughly behaves as an oscillating decaying 
exponential, with the frequency and decay rate related to the poles of the dispersion 
relation, as usual Landau damping on homogeneous backgrounds. 

ii) After this transient, the algebraic decay sets in, and the frequency changes; the 
decay exponent and the new frequency are now governed by the dominant branch point 
singularity on the real axis of the dispersion relation. 

These results prompt several questions, which this work does not answer. First, 
when the stationary distribution has a compact support in action variable, the edge 
of the support may create a singularity stronger than the bottom of the potential 
well. In some cases, one would then expect a 1/t asymptotic decay, with a frequency 
corresponding to the action at the edge of the support. However, a stationary 
distribution with compact support may also sustain purely oscillatory modes [20] , which 
do not decay at all. We have not been able to find a stationary state with compact 
support and no oscillatory mode to test the possibility of a 1/t decay. Second, the 
local extrema of the function Q(J) may also create a different type of singularity, and 
thus modify the asymptotic behavior of a perturbation. Recently, Bouchet and Morita 
have studied a similar situation in the context of the 2D Euler equation, unveiling the 
phenomenon of "vorticity depletion" close to these local extrema [30]. The possibility 
of a similar behavior for the Vlasov equation is an open question. Third, we picked 
up the slowest decaying contribution among Fourier modes with respect to the angle 
variable, and we have not investigated if the sum over the infinite number of Fourier 
modes affects the asymptotic decay. 

It is also of primary interest to understand the asymptotic behavior of a 
perturbation in a three dimensional setting. When the potential created by the 
stationary state is integrable, angle-action variables can be defined, and the method 
used in this article is viable: one would have to study the singularities of the dispersion 
relation in this case. In a situation where the potential created by the stationary state 
is not integrable, the strategy would fail. 
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Finally, our numerical study required to compute a Vlasov evolution with good 
precision for a long time, which raises difficulties. To overcome them, we have introduced 
a particle method, the weighted particles approach. It proved particularly well-suited 
for our purpose, and it would be interesting to investigate the reasons for this. 
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Appendix A. Estimate of dk(t) 

Each function afc(cu) has singularities of z u lnz type at to = mu a for all m. It may be 
written as 



where the A m (z) are supposed to be analytic, and the possible singularities of B at 
uj = mu a are weaker than z Vm In z. Isolating one term in the sum over m, we have to 
compute the following inverse Laplace transform, dropping the index m for simplicity: 



JT 

We assume that A(z) rapidly decreases for \z\ —> oo. We take a branch cut in the lower 
half of the imaginary axis as shown in Fig. IA1I 




(A.l) 



m 




I A(z — mu a )(z — mu a Y ' ha.(z — mu a )e lzt dz 



(A.2) 




Figure Al. (color online) Modification of the Bromwich contour T to the closed path 
r + + Tb + Tl in the complex z plane. The wavy line represents a branch cut 
associated to the logarithmic singularity at z = 0, which is assumed as the unique 
singularity. The numbers 37r/2 and —tt/2 are arg(z) on the present Riemann sheet. 
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To compute the integral over the Bromwich contour T, we add the paths Tr, T r 
and Tl to the Bromwich contour T, and make a closed path. From the assumption no 
singularity is enclosed by the closed path, and hence 

A(z)z"ln{z)e- lzt dz = 0. (A.3) 

r+r R +r B +r L 

Thanks to the factor e~ lzt of the integrand, contributions from Tr and Tr vanish, and 
hence 

r A{z)z v \nze- izt dz = [ A{z)z h ' \Yi{z)e~ izt dz (A.4) 



o 

+ / A(re- m/2 ){re~™ l2 y\n{re-™ /2 )e- ire ~ ml2t e- n/2 dr (A.5) 



oo 
oo 



Due to the branch cut, we have to distinguish A(re i3n ^ 2 ) from A{re l7r / 2 ) ; and we denote 
them by A^(r) and A R (r) respectively. The integral is hence written as 

POO 

A(z)z In ze~ izt dz = i{-i) u / [A L (r) - A R {r)]r u \n(r)e- rt dr 

Jo 



POO 

-if - / [3A L (r) + A R (r)y e - rt dr. (A.6) 
2 Jo 



The two functions A^(r) and A R (r) coincide in the limit r — > 0. Around the singularity 
r = 0, we can therefore estimate 

A L (r) - A R (r) oc r, (A.7) 
3A L (r) + A R (r) oc r°. (A.8) 

Using a scaling of the variable as y = rt, the integral is expressed by 

p S~1 POO 

J A(z)z In ze~ izt dz = J y v+1 {\ny - lnt) e - y dy 

r i poo 

+ ^rij o y Ve ~ Vd y- ( A - 9 ) 

The first term of the right-hand-side yields t~ v ~ 2 decay and t~ v ~ 2 In t decay, and the 
second term t~ v ~ x decay. Returning to (lA.2j) . we obtain e~ tmuJat /t u+1 as the slowest 
decay. From (jA.ip . we now see that the asymptotic decay of Ofc(t) is governed by the 
smallest u m . 



Appendix B. Weighted particles code 

The A^-body simulations are performed by the weighted particles code. In a standard 
A^-body simulations, we would prepare N initial positions and momenta by drawing 
random numbers according to a given initial distribution. An initial condition prepared 
in this way has fluctuations of order 0(l/y/~N). In the weighted particles code, we 
prepare N initial positions and momenta as lattice points of a square lattice, and give 
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to each lattice point a weight proportional to the initial distribution we want to sample. 
The concrete algorithm is as follows. 

We span the \x space with a regular lattice, having N points. The upper and lower 
boundaries ±p max for the lattice in the p-direction must be set such that the initial 
distribution f(x,p) for \p\ > p max is negligible. Let i = 1,---,N denote the lattice 
points, and (xi,Pi) their coordinates. We assign the weight Wi = Cf(xi,Pi) to each 
lattice point, where the constant C is defined by the normalization 

N 

5^ = 1. 

i=l 

We put one particle on each lattice point, and the particles move on /i space following 
by the canonical equations of motion for the HMF model 

ii(t) = Pi(t), Pi(t) = -M x (t) saxxiit) + Myit) cosxj(t), 

where the suffix i runs from 1 to N, and M x (t) and M y (t) are defined by 

N N 

M x {t) = ^Wicosxiit), M y {t) = Y,WiSm Xi {t). (B.l) 

i=i i=i 

We note that the lattice is used only to define the weight iWj, and to set the initial 
condition for the i-ih particle as (xj(0),pj(0)) = (xj,pj). It is worth stressing that 
particles are not fixed at lattice points, but move in the whole \i space, keeping their 
initially assigned weights W{. The magnetization (M x ,M y ) is computed using the 
evolving positions and the initially fixed weig ht as (IBT]1 . 

In a semi-Lagrangian code [32], the distribution is defined at fixed lattice points. To 
evolve the distribution over a time step At, one computes the inverse temporal evolution 
of a particle during At with initial condition given by the lattice (xi(0),pi(0)) = (xi,pi). 
Integrating the Vlasov equation, one gets f(xi,pt, At) = f(xi(—At),pi(—At),0). The 
point (xi(—At),Pi(—At)) does not coincide with a lattice point generally, so that the 
value of the distribution at this point is obtained by interpolation. The semi-Lagrangian 
code thus requires three discretization parameters: the time step At, and two spatial 
and velocity discretization parameters Ax, Ap. In order to insure numerical stability, 
the time step must become small as spatial and velocity discretizations become small. 
As a result, we need a small time step for computations with high spatial resolution, 
which results in heavy computations to reach the long time regime. 

In the weighted particles code, we may set the Ax and Ap discretizations 
independently of the time step At, and hence the computational burden may be reduced 
by taking a rather large At, still giving a good enough accuracy. 

The weighted particles code has further advantages against the semi-Lagrangian 
code: 

(1) The Vlasov equation has an infinite number of conserved quantities which are 
J (f(x,p)) l dxdp for / G N, but it is known that the semi-Lagrangian code cannot preserve 
them for I > 2. In the weighted particles code, the weight Wi is a fixed value and hence 
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all such quantities, approximated by ^^Wi) 1 , are preserved exactly. (2) The semi- 
Lagrangian code uses interpolation. It is a delicate step to obtain the temporal evolution 
of the distribution function, and the accuracy of the code depends on the algorithm of 
interpolation. The weighted particles code does not require any interpolation. 

Let us remark that the weighted particles algorithm is highly parallelizable and 
its convenient structure makes it possible to use a lot of parallelization methods. In 
particular, it allows to take advantage of the available computer architecture, be it 
a cluster with distributed memory or shared memory. In the case of a distributed 
memory, the mean field property allows to restrict the communication between node 
to the magnetization, which can be computed piece by piece on each node. It is thus 
possible to compute the long time evolution of the system for a very large number of 
weighted particles. 

One of the disadvantages of the weighted particles codes is that this code cannot 
compute the temporal evolution of the distribution directly. We can obtain a coarse- 
grained distribution, but its resolution is lower than the initially given lattice. Another 
disadvantage may be the limitation of objects for which the weighted particles code 
works well. The weight Wi on a lattice point corresponds to set several particles with the 
same initial condition. If these particles were given slightly different initial conditions, 
they would eventually separate as time goes by. In the weighted particles code, they 
remain together. Consequently, weighted particles code might have to be improved if it 
is to be used in order to observe more drastic changes of the distribution function, such 
as violent relaxation from a waterbag initial state to a Lynden-Bell quasi-stationary 
state. Further investigations are needed to understand why the weighted particles code 
seems to work so well in our case. The numerical tests and theoretical arguments given 
in [33] may be a first step in this direction. 



Appendix C. Extraction of a rotating part 

When we use a asymmetric perturbation such as the one given in (1661) . we observe a 
small rotation of the x and y magnetizations M x (t) and M y (t) (see Fig. IC1I) . In this 
case, it is not possible to directly define M^\t) and My~\t) by substracting the long 
time average. Our method is then to use a running average. We define M^\t) and 
M§\t) such that 

i rt+At 

M«(t) = M x (t) - — J M x (t>)dt>, (C.l) 

and 

i rt+At 

M«(t) = M y (t) - — j M y (t>)dt>. (C.2) 

In order to compute the right exponent of a power law fit, we have to choose the 
parameter At. However our different tests (see Fig. IC2p show that modifying At does 
not change much the result. We have taken At = 5 for Fig. [4j This is roughly the time 
needed to observe one oscillation of M y (t), and two oscillations of M x (t). 
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Figure CI. Slow rotation of M^p(t) (upper panel) and My X \t) (lower panel) for an 
initial condition perturbed by a sine, as in (|66[) . The numerical simulation was done 
using the weighted particles code with T = 0.1, a = 0.1 and N = 10 9 . The initial 
weighted points u>i are equally distributed in ] — w, ir] x [—3, 3]. 
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